//standard
#include <istream>
#include <numeric>
#include <vector>
#include <string>
#include <fstream>
#include <iostream>
#include <algorithm>
#include <cstdlib>

// boost
#include <boost/format.hpp>

// components
#include <tree.h>
#include <index_string.h>
#include <ultoa.h>
#include <progress_indicator.h>
#include <print_error.h>
#include <portable_timer.h>
#include <perlutil.h>
#include <sparse_matrix.h>
#include <sparse_matrix_io.h>
//#include <matrix.h>
//#include <matrix_io.h>
#include <bcpp_parsefasta.h>
#include <stlhelper.h>
#define DEBUG_LEVEL 6
#include <debug_func.h>
#include <tree_extract_sub_tree.h>
#include <tree_ops.h>

#include <boost/tuple/tuple.hpp>

#include "get_arg.h"

using boost::format;
using std::vector;
using std::ostringstream;

using std::string;
using std::fill;
using std::ofstream;
using std::cout;
using std::cin;
using std::cerr;
using std::ifstream;
using std::pair;
using std::istream;
using std::max;
using std::min;
using namespace bcpp_matrix;
using namespace bcpp_tree_ops;
using namespace bcpp_tree;
using boost::tie;

//________________________________________________________________________________________

//	index_unique_names

//________________________________________________________________________________________
void index_unique_names(const vector<string>& names, t_index_string& indexed_names)
{
	// index current names
	for (unsigned i = 0; i < names.size(); ++i)
	{
		bool added;
		unsigned index = indexed_names.index(names[i], added);
		if (!added)
			throw std::runtime_error("The matrix contained separate entries (i.e rows " +
					ultoa(index + 1) + " and " + ultoa(i + 1) + ") with identical names ["
					+ names[i] + "]");
	}
}


//________________________________________________________________________________________

//	index_unique_names

//________________________________________________________________________________________
void filter_matrix_by_leaves(	const t_sparse_matrix_f&							orig_matrix,
								const t_index_string&						orig_matrix_indexed_names,
								const std::vector<const t_tree::t_node_type*>&	leaves,
								t_sparse_matrix_f&								sub_matrix)
{

	//
	//	go through leaves
	//        add indices of matrix
	//
	vector<unsigned> leaf_indices;
	for (unsigned i = 0; i < leaves.size(); ++i)
	{
		t_index_string leaf_names;

		//
		//	make sure this valid_name has not been added before
		//
		bool added;
		leaf_names.index(leaves[i]->get_name(), added);
		if (!added)
			throw std::runtime_error("The tree contained separate leaves with identical names ["
					+ leaves[i]->get_name() + "]");

		//
		//	find leaf position in original matrix or barf
		//
		unsigned index;
		if (orig_matrix_indexed_names.is_indexed(leaves[i]->get_name(), index))
			leaf_indices.push_back(index);
		else
			throw std::runtime_error("The tree contained a leaf ["
					+ leaves[i]->get_name() + "] not found in the original matrix");
	}

	orig_matrix.subset(leaf_indices, sub_matrix);

}


void report_branches_extracted(const op_extract_branches& sub_trees, ostream& VERBOSE)
{
	// how many subtrees extracted?
	VERBOSE << "\t" << sub_trees.node_UIDs.size() << " branches extracted.\n";

	if (sub_trees.cnt_tree_extracted_by_len)
		VERBOSE << "\t" << sub_trees.cnt_tree_extracted_by_len
					<< " branches with lengths above threshold.\n";
	if (sub_trees.cnt_tree_extracted_by_size)
		VERBOSE << "\t" << sub_trees.cnt_tree_extracted_by_size
					<< " branches with number of leaves above threshold.\n";
}

bool read_tree(bcpp_tree::t_tree& tree, unsigned& cnt_lines)
{
   unsigned tree_line_pos = cnt_lines;
   cin >> t_tree_io_wrapper(tree, cnt_lines);
   if (!cin)
	   return false;
   if (!tree.get_name().length())
	   throw std::runtime_error("Tree on line " + ultoa(tree_line_pos) +
								   " has no name");
   return true;
}

void read_matrix(const string& tree_name, t_index_string& orig_matrix_indexed_names,
							t_sparse_matrix_f& orig_matrix, unsigned& cnt_lines,
							t_progress_indicator& dots, t_program_args& args)
{

    //unsigned matrix_line_pos = cnt_lines;
	//unsigned cnt_duplicate_pairs = 0;
	cin >> t_sparse_matrix_istrm_wrapper_f(orig_matrix, args.input_format,
										cnt_lines, dots);
	if (!cin)
		throw std::runtime_error("Missing matrix following tree [" + tree_name + "] "
									"on line #" + ultoa(cnt_lines));
	// steal name from tree if missing
	if (!orig_matrix.matrix_name.length())
		orig_matrix.matrix_name = tree_name;
	//if (cnt_duplicate_pairs)
	//	std_print_warning(args.VERBOSE(2), ultoa(cnt_duplicate_pairs) + " duplicate pairs"
	//						" for matrix [" + orig_matrix.matrix_name + "]"
	//						" starting on line #" + ultoa(matrix_line_pos));

	// index names
	index_unique_names(orig_matrix.names, orig_matrix_indexed_names);
}

void extract_sub_matrix(    const bcpp_tree::t_node&	node,
							const t_index_string& 			orig_matrix_indexed_names,
							const t_sparse_matrix_f& 				orig_matrix,
							t_sparse_matrix_f&					sub_matrix,
							const string&					sub_matrix_name)
{
	// gather leaves
	bcpp_tree_ops::op_extract_branch_leaves tree_branch(node);

	// get sub matrix
	filter_matrix_by_leaves(orig_matrix, orig_matrix_indexed_names,
								tree_branch.leaves, sub_matrix);
	// submatrix name
	sub_matrix.matrix_name = sub_matrix_name;

}

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
int main (int argc, char *argv[])
{
	debug_func_cerr(5);

	// we use cout to print progress so flush per call
	//cout.setf(std::ios_base::unitbuf);

	std::ios::sync_with_stdio(false);

	t_program_args args;
	try
	{
		if (!process_command_line_options(argc, argv, args))
			return 1;


	//    char buff[4096];
	//    cin.rdbuf()->pubsetbuf(buff, sizeof buff ); // enable buffering

		ofstream unchanged_tree_file(args.unchanged_tree_file.c_str());
		if (!unchanged_tree_file)
			throw std::runtime_error("Could not open the file \"" +
										args.unchanged_tree_file + "\"");

		t_progress_indicator dots(args.VERBOSE_ONLY(1),  1000);
		unsigned cnt_trees = 0;
		unsigned cnt_lines = 0;
		args.VERBOSE(1) << "\tExtracting large tree branches...\n";
		unsigned sub_tree_cnt = 0;
		while (cin)
		{
			//
			//	read tree from cin
			//
			args.VERBOSE(2) << "\n\tReading tree...\n";
			bcpp_tree::t_tree tree;
			if (!read_tree(tree, cnt_lines))
			{
				args.VERBOSE(2) << "\tNo more trees.\n\n";
				break;
			}

			//
			//	Read corresponding matrix
			//
			args.VERBOSE(2) << "\n\tReading matrix...\n";
			t_sparse_matrix_f							orig_matrix;
			t_index_string						orig_matrix_indexed_names;
			read_matrix(tree.get_name(), orig_matrix_indexed_names, orig_matrix,
									cnt_lines, dots, args);

			++cnt_trees;

			//
			//	extract heads of branches
			//
			args.VERBOSE(2) << "\tSelect large tree branches above thresholds...\n";
			op_extract_branches sub_trees(	tree, args.node_limit, args.len_limit);
			if (!sub_trees.node_UIDs.size())
			{
				args.VERBOSE(2) << "\tTree unchanged...\n";
				unchanged_tree_file << tree << "\n";
				continue;
			}
			report_branches_extracted(sub_trees,  args.VERBOSE(2));



			//
			//	write denuded tree in nh format and extract sub-trees corresponding to
			//			truncated branches
			//
			//		add offset when referring to sub-branches so that different trees
			//			on the same file have unique branch names
			if (!tree.get_name().empty())
				cout << ">" << tree.get_name() << "\n";
			op_nh_strm_post_branch_extract main_tree_after_extract(	cout,
																	sub_trees.node_UIDs,
																	sub_tree_cnt);
			main_tree_after_extract(tree.get_root());
			if (main_tree_after_extract.cnt_sub_trees != sub_trees.node_UIDs.size())
			{
				cerr << "help\n";
			}

			//
			//	Extracting branches
			//
			args.VERBOSE(2) << "\tExtracting branches...\n";
			for (unsigned i = 0; i < main_tree_after_extract.cnt_sub_trees; ++i)
			{
				args.open_output_files(sub_tree_cnt + i + 1);

				//
				// print sub-tree
				//
				if (args.f_tree_file)
				{
					// tree name
					args.f_tree_file
						<< ">" << main_tree_after_extract.sub_tree_name(i) <<  ' '
								<< args.description << "\n"
						<< *main_tree_after_extract.nodes[i]
						<< "\n";
				}


				//
				// print sub-matrix
				//
				if (args.f_matrix_file)
				{
					// get sub matrix
					t_sparse_matrix_f sub_matrix;
					extract_sub_matrix(*main_tree_after_extract.nodes[i],
										orig_matrix_indexed_names, orig_matrix,
										sub_matrix,
										main_tree_after_extract.sub_tree_name(i) + ' ' +
											 args.description);

					args.f_matrix_file << t_sparse_matrix_ostrm_wrapper_f(sub_matrix,
																	args.output_format,
																	dots);
				}
			}
			sub_tree_cnt += main_tree_after_extract.cnt_sub_trees;

		}
		if (!cnt_trees)
			throw std::logic_error("No tree in New Hampshire format was found.");
	}
	catch(std::exception& e)
	{
		std_print_error(args.CERR(), e.what());
		return 1;
	}

	args.VERBOSE(1) << "\tAll finished!!\n\n";

	return 0;
}



